package scu.maqiang.piezoelectricity;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L2PE;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class Piezoelectricity1DL2 {

	public static void main(String[] args) {
		Mesh1L2 mesh = new Mesh1L2();
		mesh.line(0, 1, 128);
		
		FES1L2PE fs = new FES1L2PE(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		
		double G1 = 1.0, ep1 = 0.5, e1 = 0.1;
		fs.assembleStiff(new double[] {G1, ep1, e1}, BVPType.COMMON, A);
		double f1 = 1.0, rhof = 2.0;
		fs.assembleSource(new double[] {f1, rhof}, BVPType.COMMON, RHS);
		double[] x = new double[fs.GetNdof()];
		fs.applyBC_MBN(A, Direct.All, 1);
		fs.applyBC_MBN(RHS, Direct.All, 0.0, 1);
		//A.sort();
		//System.out.println(A);
		//System.out.println(MVO.toString(RHS));
		
		DirectSSolver solver = new DirectSSolver(A);
		x = solver.LUSolve(RHS, 0, 1e-6);
		//NewIterSSolver solver = new NewIterSSolver(A);
		//solver.GMRES_SSOR(RHS, x, 1.5, 20, 1);
		//A.save("matFile.txt");
		//MVO.saveArray("RHSFile.txt", RHS);
		//solver.CGNE(RHS, x, 1);
		//Arrays.fill(x, 0.0);
		//solver.CGNE(RHS, x, 1);
		System.out.println("x min:" + MVO.min(x) + "max: " + MVO.max(x));
		//Arrays.fill(x, 0.0);
		//solver.CGNR(RHS, x, 1);
		//System.out.println("x min:" + MVO.min(x) + "max: " + MVO.max(x));
		double[][] u = new double[4][mesh.getNv()];
		fs.extractComponent(x, u);
		
//		//ѹ�����⾫ȷ��
		ScalarFunc uExact = (xy, label, param) -> 
		                        (rhof * e1 - f1 * ep1)/(G1 * ep1 + e1 * e1) * (0.5 * xy[0] * xy[0] - xy[0]);
		ScalarFunc phiExact = (xy, label, param) -> 
		                        -(f1 * e1 + rhof * G1)/(G1 * ep1 + e1 * e1)* (0.5 * xy[0] * xy[0] - xy[0]);
		
		u[2] = fs.valueFromFunc(uExact, null);
		u[3] = fs.valueFromFunc(phiExact, null);
		mesh.toTecplot("PESolution.dat", u);
	}

}
